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\ Abstract 

F ^ | A flexible spectral code for the study of general relativistic magne- 

tohydrodynamics is presented. Aiming at investigating the physics of 
slowly rotating magnetized compact stars, this new code makes use of 
various physically motivated approximations. Among them, the relativis- 
tic anelastic approximation is a key ingredient of the current version of the 
code. In this article, we mainly outline the method, putting emphasis on 
algorithmic techniques that enable to benefit as much as possible of the 
■ non-dissipative character of spectral methods, showing also a potential 

' astrophysical application and providing a few illustrative tests. 

o ' 

1 Introduction and motivation 

f^*) . During the last decades, Astrophysics has seen fast progress in its ob- 

servational techniques, and a very wide flow of data is now coming from 
^-l several orbital and ground-based detectors. This huge amount of infor- 

mation is a challenge not only to the data analysis community, but also 
to theorists, as accurate numerical models become more and more re- 
quired to explain various astrophysical observations and to put to the 
test more qualitative theoretical predictions. Among the processes that 
k>( ' are the most interesting and complicated from both physical and nu- 

merical points of view, those related with strong gravitational and/or 
magnetic fields play crucial roles inside or in the vicinity of compact 
objects. Hence, they are involved in the description of numerous phe- 
nomena, such as are massive stars core collapses, emergence of jets from 
active galaxies, emission of gamma-ray bursts (GRB) or of gravitational 
waves (GW), merger of binary neutron stars (NS) or just evolution of 
isolated NS. As a consequence, various codes for r elativistic Magneto - 
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HydroDynamics (MHD) ha v e been developed [e.g . 
De Villiers fc Hawlevl (120031). iGammie et al\ d200j>. 



Koide et al 



Komissarov 



1999) 



2004) 



Shibata fc Sekiguchil (|2005l ). iDuez et al\ (|2005l ), lAnt6n et ql.l |2006l ) J^frzuno et all 
(|2006l )]. with different scientific purposes, and therefore distinct choices 
in the physical approximations, in the writing of the equations and in the 
numerical implementation. 

In the following, we present a new relativistic MHD code, mainly de- 
scribing the methods and some basic tests. This code is tridimensional 
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and based on spectral met hods, quite similar to the r elativistic hydrody- 
namical code introduced in Villain & Bonazzola (2002) (hereafter referred 
to as Paper I) and IVillain. Bonazzola fc Haensell (|2005h (Paper II). The 
new code was designed to be flexible and to fulfill several demands, the 
first of them being to solve MHD equations in the fixed curved space- 
time associated with the steady-state configuration of a slowly rotating 
NS (i.e. it currently works in the Cowling approximation). Those equa- 
tions are solved in the whole interior of the star (including the origin), 
using spherical vector components and spherical coordinates, which is al- 
most necessary to easily implement various boundary conditions (BC), 
and hence to be able to deal with a wide spectrum of physical situations. 
Another requirement for exploring numerous physical phenomena is that 
of the code being robust and stable in the sense of allowing to follow the 
evolutionary tracks during many rotational periods. In such a way, it 
can also be used to study the stability of equilibrium configurations of 
rotating magnetized NSs. As we shall illustrate, this is made possible by 
the non-dissipative character of spectral methods and by some cautions 
taken in the algorithm, such as having the numerical magnetic field that 
remains divergence- free. 

The article is divided in several sections: in Sect. [2] we summarized 
the typical orders of magnitude and timescales involved in the physical 
situations that we shall study in the future, using the resulting approxi- 
mations to write the relativistic equations of motion in a Newtonian-like 
way. Then, the numerical method is described in Sect. [3] We discuss some 
possible applications in Sect. [4] Finally, Sect. [5] contains tests performed, 
whereas final conclusions are gathered in Sect. [6] 



2 Timescales and equations of motion 

We shall briefly remind orders of magnitude and timescales for the type 
of problems to be addressed with the code: We are interested in the MHD 
of old but not too cold (not superfluid) compact stars, which have radii 
of some 10 6 cm, mean baryonic densities n ~ 10 14 g/cm 3 , and magnetic 
fields of 10 12 to 10 16 gauss. We consider only compact stars whose density 
does not vanish at the surface (the so-called bare strange stars) or which 
have a rigid crust. In the latter case the domain of integration is limited 
to the liquid part of the star, completely liquid stars (just born very hot 
NSs) being much more difficult to treat since BC of wind solutions have 
to be imposecQ. 

Most of known isolated rotating NSs do not have rotational frequency 
larger than 300 Hz (~ fi = 1880 rad/s), which gives a period T rot = 
3.3 10 -3 s. Retaining this typical value implies that to keep only terms 
linear in Q is a good approximation, be it for the relativistic Euler equation 
or the background metric. Using also the conformal flatness approxima- 
tion, the latter is written as in Papers I and II in a (3 + l)-like isotropic 
way (in c = 1 units) 

ds 2 = - N 2 dt 2 - 2r 2 sin(60 2 a 2 N 3 d(j>dt + a 2 dl 2 , (1) 

1 This kind of problem can be treated with two different numerical techniq ues: A finite 
differe nce method [like the one developed by the Valencia group, see for instance Anton et al. 
1(20061 )] during the in-falling an d cooling phase, then the present algorithm once the crust is 
formed. On NSs winds, see also Bucciantini et al. (2006). 
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where dl 2 is the Euclidean infinitesimal 3-length, while r and 8 are spher- 
ical coordinates, with the lapse N, the conformal factor a, and the shift 
N3 that only depend on r. 

Another important velocity for magnetized NSs is the Alfven velocity 
defined as c a ~ B/^/Tlm ~ 10 4 B 12 cm/s [with B X2 = B/(10 12 gauss)], 
leading to the Alfven crossing time tai = 10 2 B-jT 1 s (and to Vax = 
10~ 2 B12 Hz). This time and the period are much larger than the acous- 
tic crossing time r ac , whose value comes from the sound velocity that, 
in compact stars, is approximately equal to one third of the light ve- 
locity (c s ~ 0.3 c — * T ac ~ 3.3 10~ 5 s, with corresponding frequency: 
Vac ~ 3. 10 5 Hz). As Tac is much shorter than the other characteris- 
tic times, we shall filter acoustic waves through the so-called relativistic 
anelastic approximation introduced in Paper I: 

4=ft(V57' i ) = 0, (2) 

with e the determinant of the Euclidean 3-metric and V 1 = nNa 3 !! 1 , 
where U l is the spatial part of the fluid 4- velocity (i £ [1,3]). Notice 
that contrary to what was done in Papers I and II, this relation is not 
necessarily linearized here, for we shall also use the divergence-free char- 
acter of the P % 3-vector as a key feature for the algorithmic (see Sect. 
[3]). Moreover, we shall already mention that the momentum density V 
(which is the 1-form associated to V) is the main variable retained in the 
hydrodynamical part of the new code, while in previous works (Papers 
I and II), this variable was the velocity V. The reason of the previous 
choice was that the vanishing matter density at the surface of the star (at 
least in Paper I), would not have made possible in general to recover the 
velocity from the momentum densitjQ. 

As the star can b e described as composed of a single neutral fluid (see 
Paper II), we follow IWilsonl |l972h and use both the form of the metric 
|T]) and the anelastic approximation ([2]), to write Euler equation in the 
Newtonian-like way 



dtVi + 



f J \n V / 
VuVv =» , a 3 N 2 n 



Viff^ + J» , (3) 



2a? NT aM / 
where / = p + P, whereas Ti^ is the Lorentz force. 

Mainly due to the frame-dragging effect, to curvature of space-time 
and to the possible non-ideality of the plasma, the Lorentz force contains 
numerous terms that we shall not all make explicit here. We shall just 
write it in such a way to make more visible the equivalent of the Newtonian 
force, which will be used in Sect. [5] to proceed to basic tests of the code 

Fiy, J 1 " ~ (V A B) A B + RT, (4) 

47T x ' 

where Z is a metric factor, B a magnetic field-like vector and RT en- 
compasses all other Relativistic Terms (related with frame-dragging, non- 



2 The division of V by the matter density n is numerically possible only for stiff equation 
of state (r > 2). 
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ideality of the conductivity, etc.). B is defined as 

B> = \ri 3k Fj k , (5) 

with rfi h the tridimensional Euclidean Levi-Civita volume form and Fjk 
the Maxwell electromagnetic tensor. Its evolution is ruled by the so-called 
induction equation, which comes from relativistic Ohm's law[f| applied to 
a neutral plasma 

kJ„ = F^U" , (6) 

where k is the conductivity^, and from the homogeneous Maxwell equa- 
tions 

d a F M + dp F ja + d 1 F afj = 0. (7) 

The relativistic induction equation for B is written in a Newtonian-like 

way as 

d t B + VA (BA V) + VA (nj) = 0, (8) 

where V 1 = U z /U° and J % = N ,P . Notice that due to the frame-dragging 
effect, ,P is not directly related to the curl of B by the inhomogeneous 
Maxwell equation, but we shall not enter more into the detail here, and 
we just mention that while testing the code, we used the Newtonian ex- 
pression of Faraday's equation in the MHD approximation 

i J = V A B. (9) 

47T 

The main reason behind that approximation is that it does not change 
the algorithmic but allows to verify the conservation of energy (see Sect [5} 
and consequently the whole mechanics of the code and its non-dissipative 
character. Furthermore, we would like to insist on the fact that what we 
present here is the simplest version of the code, containing the minimal 
number of physical assumptions. Thanks to the required flexibility, more 
physical situations can be easily considered (for example strong stratifi- 
cation, see further), by following the same strategy as going from Paper 
I to Paper II. However, notice that when treating non-linear problems, 
we will need some kind of artificial viscosity in order to smooth the so- 
lutions. Yet, the main difference between spectral and finite difference 
scheme codes is that with spectral codes, artificial viscosity and resistiv- 
ity are put as differential operators, implying that we perfectly know their 
behaviour, while they also help to impose boundary conditions. Since for 
linearizeqj problems, we do not need, in principle, dissipative terms, spe- 
cific problems can arise due to the fact that instabilities have much more 
possibilities to develop. Consequently, the linearized non-dissipative ver- 
sion of the code is still in progress, even if some tests are already provided 
here (see Sect. 0J, showing that the global strategy works well. 



3 See for instance iLic hncrowicz J 1967f) ■ 

4 Notice that the decay time of magnetic field in NS is, as the viscous time, much larger 
than t he dynamical time under consideration here. See for instance Goldreich & RciscncEgcr 
l|l992f ). 

5 Notice that what we call the linear case is a linearization around the rigid rotation solution, 
implying equations of motion very similar to the ones described in Paper I and II. 
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3 The method 



3.1 Implementation of the anelastic approxima- 
tion 

As in Papers I and II, the solving of the equations (linear or not) relies 
on two main ingredients: 

• the Helmholtz theorem, that is to say the decomposition of vectors 
into divergence-free components and potential parts, 

• the use of so-called angular potentials. 

However, working with as variable the momentum makes the imple- 
mentation of the anelastic approximation quite different, and we proceed 
by considering the equation of motion (Eq. [3| written in the following 
way: 

dtPi + d t P = Si , (10) 

where P is a kind of geometrically modified pressure (P = a TV 2 II with 
dll = ndP/f), while the term Si represents all other terms appearing 
in Eq. Q (non-linear included) and which are considered as a known 
source. Applied to this known source term Si, Helmholtz theorem gives a 
divergence- free component Si and a potential part di$f: 

Si = Si + diV. (11) 

Note that the above decomposition is not unique, since one can add 
an arbitrary harmonic function to the potential $ (or equivalently 
subtract its gradient to Si) without altering the divergence of Si. As we 
shall see in the following, this degree of freedom is used to impose an arbi- 
trary BC for V (for example P r = 0). Due to the anelastic approximation 
and the exact compensation that is required between diP and di*S>, the 
equation to be solved reads 

d t Vi = Si. (12) 

If the non-linear equations are considered, a viscosity £ (supposed to 
be constant) is added, and we have 

dtVi + CAVi = S i (13) 

Once P is obtained (P = ^f), the new source term for the next time 
step can be calculated (in the barotropic case) by using the equation of 
state and the frozen metric. But before entering more into the detail of 
solving such equations (the induction equation being technically identi- 
cal), we shall now briefly recap on angular potentials and their use in 
spherical geometry. 

3.2 Vector angular potentials 

As mentioned in Sect.[T] we are interested in solving vectorial partial dif- 
ferential equations (PDE) employing spherical components and spherical 
coordinates. The advantage of this choice is the possibility to impose 
exact boundary conditions at the surface of the star. The drawback, of 
course, are the coordinate singularities appearing in the PDEs, but with 
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spectral methods those difficulties can be ea sily overcome if we intro- 
duce the angular pote ntials r/ and [i defined as l|Bonazzola fc Marcklll990l . 
iBonazzola et al\\l99di ) 



Be = der/ —rd<t,fj,, B$ = -^—^d^r) + dg^ (14) 

and smO 



or analogously 

d 2 r] cosddr) 1 . , 

A e ,^=^ + — - + — (15) 

_dBe_ + cos6_ B _]_dB ± 
88 sin # 8 sin # 90 
Now, if we expand r\ in spherical harmonics we have 

A 9 , 0?7 = -1(1 + 1)1] (16) 

Analogously for /j 

A e , 0A i = <9 9 B + - -J— 9 B 9 (17) 

To illustrate the use of these potential by some examples, we shall 
mention that the divergence of a vector B can be written 

divB = d r B r + -B r + - (dgBe + ^rB e + -J—d^) , (18) 
r r \ smO sm& / 

which reads, in spherical harmonic representation, 

d r B r + -B r - -1(1 + l)r) . (19) 

r r 

In the same way, to deal with an equation like AB = S for a divergence- 
free vector B, the decomposition in potentials of both B and S leads to 
solve ^ 

" " L!r + ~ir B - +\(2-l(l+ 1)) B r = Sr , (20) 



dr 2 r dr 
d 2 ri 2 dr? 1 
dr 2 r dr r 



2+-7Z + Zs(-Kl + ^)v + 2B r )) = Vs , (21) 



d 2 a 2 da 1,„ „ . 

where Eq. (|20|l was obtained using Eq. (|19p and the divergence- free nature 
of B. Notice that the equation for the toroidal part fi is decoupled from 
the other two equations. 



Having obtained ordinary differential equations for the radial variable, 
we see that the problem of the singularities on the axis 6 — ~tt, 8 = n is 
easily solved by means of spectral decomposition on spherical harmonics. 
Similarly, problems of singularities at r = are solved by expanding 
the solution on a complete set of functions that have a good analytical 
behaviour at r — (e.g. Chebyshev polynomials or a linear combination 
of them, see Appendix of Paper I) . 
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3.3 Divergence-free decomposition 

We shall now describe in more detail the new method that is used to apply 
Helmholtz theorem and decompose a vector field into a divergence- free and 
a potential parts. Consider again the equations of motion Eq. (|12l I13p . 
where Si is a source term that is supposed to be known at the time 
tj-i, and that we have to decompose into a divergence- free part Si and a 
potential part di&, 

Si = Si + ft$ . (23) 

The ordinary way to make this decomposition is described in all analysis 
textbooks and was applied in Paper I: one simply takes the divergence of 
the Eq. f|23p and solves the Poisson equation 

A$ = divS. (24) 

Once <!> is obtainecjf], we got Si by Si — Si — di^ . 

The problem in the numerical implementation of this method is that 
in order to obtain Si from <!> one needs to compute its gradient, while the 
numerical computation of a derivative means higher round-off errors and 
generation of numerical instabilities. To avoid this problem, in the MHD 
code, we directly get S by first computing C v , the r\ angular potential of 
the curl of S, and then solving the system 

(curlS) v = {curlS) v , divS = . (25) 
After an expansion in spherical harmonics, this system reads 

dr, 77 S r r dS r 2 - 1(1 + 1) n (0(f . 

— H =C V , —. h -S r V = 0, (26) 

dr r r dr r r 

where rj is the potential of S. An expansion in Chebyshev polynomials 
leads to an algebraic linear system whose matrix is a 2(N r — 1) x (2iV r — 1) 
one (iV r being the number of coefficients in the expansion) that can be 
reduced easily to a 8-diagonals matrix, and in this way, S is solved without 
using $. 



3.4 Solution for the equation of a divergence-free 
vector 

In the general case, the MHD equations are not linear and the solutions 
with zero viscosity (infinite Reynolds number) or/and zero conductivity 
(infinite magnetic Reynolds number) present discontinuity in the deriva- 
tives (infinite shear). Consequently, they cannot be described by numeri- 
cal solutions with finite resolution and require specific techniques. Hence, 
in the following, we mainly consider the general case with both conduc- 
tivity a and viscosity ( not equal to zero, particular cases for which £ 
or/and a vanish (i.e. linearized problems) being only partially discussed 
in the next Section, but mainly kept for another article. Furthermore, we 
describe here only the method used to solve the equation for the magnetic 
field B, the case of the equation of the density momentum V being simi- 
lar since the dispersive equation for the magnetic field B [Eq. ©] can be 
written 

d t B + Airccurl{o curlB) = S , (27) 

6 As already mentioned, <I> is defined within a harmonic functions <3>yi that is kept as an 
additional degree of freedom used for the implementation of BC. 
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where S contains all other terms. 



To proceed, we use the decomposition of the vector B in radial com- 
ponent B r and two angular potential r\ and fi as explained above (Eg, 1141) . 
The source term is calculated at time tj+1/2 using the previous values at 
time tj and tj-i by 

sj+1/2 = _ (2g) 

For the moment, we will consider a constant, the non-constant case 
being technically identical (see Appendix of Paper I). Using the identity 

curl curl B = — A B + grad [divBj 

and taking into account the divergence-free nature of B, an expansion in 
spherical harmonics gives the system [see Eqs. ()20l21l22j )] 



d t B r + a 
dtj] + a 



+ ^(2-1(1 + 1))* 
r dr r2 



dr' 2 



d 2 r, + 2dri 
dr 2 r dr 



dtfJ, + o 



d u 2 dfi 



(29) 

(30) 
(31) 



dr 2 r dr r 

where S v and s M are respectively the angular potentials of the source 
Remember that the condition divB = reads 



^ + ^B r --l{l + 1)77 = 0. 

dr r r 



(32) 



Notice also that the equations for B r and fi are decoupled and behave like 
the Poisson equation. 



There are different strategies to solve the above system: we can, for ex- 
ample solve the Eq. (|29[) and compute 77 by using the divergence [Eq. (|32p] : 

1 / dB r „_ \ 

This method is very simple, but its drawback is that we have to perform 
a derivative with respect to r of B r , that generates high spacial frequency 
noise, which in turn produces numerical instabilities in a MHD code. We 
can also use the method described in Paper I, but we present here a dif- 
ferent one that is more robust. 



This new method is the following: After the time discretisation, a sec- 
ond order (in time) implicit scheme consists in re- writing the Eqs. (|29II30[) 
in the following way: 

~d 2 B^ 

: . , 

r 2 



3+1 



St 



di 



rl ArtW + 1 1 

+ -^- + -(2-!(l + l))B?^ 



dr 



B ] r + 



^ + ^ + 4(2-^ + D)^ 
dr 2 r dr r 2 



+ StS r 



j+1/2 



(33) 



J+ 1 



St 



rf a 



, St 
rf + -a 



d 2 rf +1 + 2 drr> +1 
dr 2 r dr 



d 2 rf 2 drf 1 . , 

—rj- H 4- + —(2B r 

dr 2 r dr r 2 



+ l(2Bl +1 - 1(1 + l) n 

Ki + i)v j ) 



+ StS; 



(34) 
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where St is the time increment and S^ +1 ^ 2 is the value of the source 
extrapolated at time tj +St/2 as explained by Eq. (|28p . Note the presence 
of the singular term (2B 3 r - 1(1 + \)rf)jr 2 at the RHS of the Eq. p4|) . 
An elementary study of the analytical properties of the above system of 
equations shows that the solutions vanish at least as r 1 " 1 at the origin 
Consequently for I > 3, the solutions vanish fast enough to avoid 
actual singularities in the apparently singular terms of Eq. (|34p RHS. For 
I = 1,2, those terms are indeed singular, but the solution is not singular 
because they do compensate exactly. However, in order to overcome this 
numerical difficulty, we re-write Eq. (I34p by using Eq. (|32p as follows: 



J+i 



St 



J] G 



, St 

rf + -a 



dV +1 

dr 2 



2 drf 
r 



3+1 



dr 



^(Ki-iW +1 ) 



d 2 r r dr r 2 



+ (S 1 + S' v 



St, 



(35) 



where S 1 is a new source term 
2 



S 1 



a 
2^ 



l + l 



dBj +1/2 

dr 



+ 2B 



3 + 1/2 



+ 2B 



■3 + 1/2 



It easy to see that for I = 1 or I — 2, since B l r = bir 1 1 + 
rf = eir 1 " 1 + 



(36) 
and 



all the terms in Eq. (|35p are regular. 



Then, the way to proceed consists in first, getting homogeneous solu- 
tions by doing the following: 

a) - Compute an homogeneous solution of the system of Eq. (|33l 134ft : 
solve with S 3 r +1 ^ 2 = for calculating B} +1 , let be this solution. 

b) - Compute an homogeneous solution r\h of the Eq. ((35jl : solve with 
S 1 = and S J V +1/2 = 

c) - Compute a particular solution rj hl of the Eq. (|35p using 5 ,J + 1 / 2 = 
and S 1 , the latter being computed from B^ 

d) - Impose, by using the homogeneous solution r/h, that the divergence 
given by the Eq. (|32ft vanishes at the boundary (r — R). 

In such a way, we have a homogeneous solution of the system (|33l 
I34p that satisfies exactly the condition of vanishing divergence (St +1 ^ 2 = 
0, 5^ +1//2 = 0) at the boundary. All those homogeneous solutions are 
stored and have to be computed at each time-step whether 5t a depends 
on time. The complete solution is computed by obtaining, in a similar way, 
a particular solution and by using the homogeneous solution to satisfy the 
BC. The condition divB — 0, as was already said, is fulfilled exactly by 
construction at the boundary and approximately in the interval < r < 
R. For a given number of degrees of freedom N r in r, the error e depends 
on the value of St a. In practical cases, we had 1CT 11 < e < 10" 5 . 



4 Non-dissipative MHD 

In the previous section, we have discussed the general MHD problem 
(non- vanishing viscosity and resistivity). As explained, in this case, the 

7 The value I = is forbidden, because we treat a pure vector field. 
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presence of viscosity and resistivity allows us to easily impose BCs and 
help to stabilize. Yet, there are astrophysical steady-state problems for 
which we really need a non-dissipative code, for example while probing the 
stability of hydromagnetic configurations. For such a class of problems, 
linearized MHD has to be used, viscosity and resistivity being no more 
necessary to smooth the solutions, being on the contrary to be avoided. 
Hence, the implementation of BCs has to be done in a different way. This 
technology is in progress, and as a first illustration, we present here an 
example of application inspired by a model of 7-ray bursts proposed re- 
ce ntlv by iPaczvriski fc Haensell (I2005T) and based on previous proposals 
bv lKluzniak fc Rudermanl (|l998h and lDai fc Lul (|l998l ). 

This model consists in a just born differentially rotating and magne- 
tized NS, that during the cooling time (and consequently the shrinking 
time) undergoes a phase transition and becomes a differentially rotating 
stratified quark star. At this point, the poloidal magnetic field is stretched 
by the differential rotation and a fraction of the rotational kinetic energy 
related with the differential rotation is transformed into magnetic energy 
of a toroidal field. Orders of magnitude estimation shows that for a star 
initially rotating at 200 Hz, with an amount of some few percent of the 
total rotational kinetic energy due to the differential rotation, energy bud- 
get is met. Moreover, a strong toroidal magnetic field (~ 10 gauss) is 
generated. If there are no losses, the magnetic field reaches a maximum 
and starts to decrease, transforming its magnetic energy back into differ- 
ential kinetic energy (at constant kinetic momentum). The duration of 
this twisting period depends on the intensity of the preexisting poloidal 
magnetic counterpart B p (~ 100 rotation periods for B p = 10 13 gauss). 
From a numerical point of view, the winding time of the toroidal magnetic 
field Bphi , takes about hundred rotation periods of the star, meaning that 
it is crucial to avoid important energy losses. 

The low dissipation condition cannot be fulfilled with the code de- 
scribed above reasonable numerical viscosity and resistivity. Hence, we 
modified it, taking into account the particularity of the problem in order 
to eliminate numerical dissipation. The strong stratification of the star 
suggests the approximation V r = and Vg = 0, since motions of the mat- 
ter inside the star tend to be parallel to isopotential surfaces: V r ~ 0, 
Vg ~ 0. Consequently, only evolutions of the toroidal components of the 
magnetic field and of the velocity V4, are pertinent to our problem. 
Moreover, as a first step, we consider in the following an incompressible 
fluid and Newtonian gravity. 

4.1 Evolution of the toroidal component 

The equations of motion for the toy model are 

d t B = curl (V A B) , (37) 

and 

d t V = — (curLB) A B . (38) 

Taking into account the assumptions described above, the axisymme- 
try and the equatorial symmetry, the equations to be solved read 

dtB^ - B r (drV^ ~ -vA ~ — (dgV^ - -^vA = , (39) 
V r J r \ r smd J 
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v 4-jrn 



(d rJ 



B r ( d r Bt + -bA + ?l (deB^ + ^B^ 
r J r \ smO 



= 0, (40) 



from which we see that in the case of rigid rotation, V$ = fir sin#, B^ 
does not evolve. 

In what follows, we shall assume in addition that the poloidal magnetic 
field can be written (taking into account the symmetries of the problem) 

B r = cos 9 (61 + b 2 r 2 cos(26») + ...), 

Be being obtained by the condition div_B = 0. 

We have to impose BCs in order to have an unique solution. A study of 
Eqs. (|39|l and (I40p shows that this system is hyperbolic, and consequently, 
one BC can be imposed at the surface of the star. Physically, what occurs 
is that a radial Alfven wave is created in order to satisfy the given BC. 
Since in the local approximation, the phase velocity Vai of an Alfven wave 
verifies 

v 2 - w2 - (k'-^) 2 41 

Al k 2 47m k 2 ' 
for an Alfven wave propagating radially, the phase velocity is 

1 -B 2 , (42) 



47rn 

which is maximum at the pole (0 — 0) and vanishes for 6 — n/2 (the equa- 
tor), meaning some possible numerical troubles due to the large variation 
of the timescale. 



Hence, the technique to solve this system is to work in configuration 
space in 9 and in the Chebyshev space for r. Once the discretisation of 8 
is done, we are confronted with Ng hyperbolic systems that can be solved 
implicitly. Since the Alfven velocity can be very small for 6 close to tt/2, 
the number of degrees of freedom N r for r must be high enough to de- 
scribe the Alfven wave that propagates very slowly along the radius. In 
conclusion, N r depends on Ng. Note that for 6 — tt/2 (the equator) the 
system is completely degenerated, and B^ vanishes as it is anti-symmetric 
with respect to the equator. 

To conclude, we emphasize that the finale solution of the system 
strongly depends on the chosen BCs. For example, if we choose B^ = 
j r =i?, the Poynting vector vanishes at the surface, and the energy is 
conservecjf]. Thus, in the next Section, we shall provide some examples 
with conditions that allows us to easily check the energy conservation and 
consequently to test the accuracy of the code. More precisely, we shall 
use for the non-dissipative case the condition B r — \ r =R at t = 0, which 
makes the operator degenerated at the surface, implying that no BC is 
required but the Poynting vector does vanish at the surface. 

Note that the radial component J r of the current vanishes at the surface too: 

1 , cos 6 

j r = -(a e s^ + — — . 

r smo 
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5 Tests 

In the following, we shall illustrate the stability of the code, but mainly 
with tests that concern linear (or linearised) problems. Yet, we would like 
to point out again that linear problems, when using spectral methods, are 
the most severe ones. Indeed, for a linear problem, a spectral code has to 
work with vanishing viscosity and conductivity, and since spectral meth- 
ods have no intrinsic numerical dissipation, numerical instabilities have 
greater opportunity to grow than in the non-linear case, where viscosity 
and resistivity can be increased until the code becomes stable. This is 
the reason why the stability conditions illustrated in the following were 
obtained in a non-trivial way, that will nevertheless be explained in detail 
somewhere else. 



5.1 Dissipative case 

Since the dissipative case is far being the most interesting one from the 
numerical point of view, we only provide one example of evolution in the 
case of a hydrodynamical coupling between an r-mode (lower frequency) 
and / = 2 f-mode (higher frequency). The calculation pictured on Fig[T] 
shows the time evolution of the component Vcf, of the velocity on the sur- 
face of the star. The calculation is done with the number of coefficient in 
spectral expansion equal to 24, 6 and 4 in r, 8 and <j> directions, respec- 
tively, and the time-step is 10 -2 (in the computational units in which the 
period of rotation is equal unity). The dimensionless viscosity parameter 
equals 0.001. 



0.004 

y 0.002 
| 

i 

> 

-0.002 



100 200 300 

no. of rotation periods 

Figure 1: Time evolution of the velocity v<f, of a given point on the surface of 
the star for a dissipative code with excited I — 2 r- and f-modes. 




5.2 Non-dissipative case 

We also present a couple of illustrative tests of the spectral code, that 
were performed with the linear version, neglecting all dissipative terms. 
On Figs. [5] and [3] the case of a pure hydrodynamical coupling between 
an r-mode and / = 2 f-mode is shown. As in previous subsection, the 
number of coefficient in spectral expansion is again 24, 6 and 4 in r, 9 and 
(j> directions, respectively, while the time-step is 10 -2 . 
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Figure 2: component of the velocity for a given point on the surface of the 
star with excited I — 2 r- and f-modes, versus time. 




500 

no. of rotation periods 



Figure 3: The same as on Fig. but for a longer run, in order to illustrate the 
lack of dissipation and the long-term stability of the code. 



Last but not least, we shall also plot the conservation of energy in 
a linear MHD code using the no-radial-magnetic field condition (see end 
of previous Section) and no dissipation term. The chosen model is ax- 
isymmetric (24 coefficients in r-direction and 16 in ^-direction) with the 
time-step equal 10 -3 and with an initial profile r sin# — r 3 sin 3 9 and an 
initial magnetic field with only a radial component: B r — smd(l—(r/R) 2 ). 
On Fig. it can be seen that for 2.5 rotation periods (2500 time-steps), 
energy is conserved up to the value of 10 -8 (the difference between the 
total initial and final energy), meaning a relative variation around 10 -5 
only. As expected, at the end of the cycle, the star has gained strong 
differential rotation and a poloidal component of the magnetic field. 



6 Conclusions 

We have described a code based on spectral methods able to handle 3D rel- 
ativists MHD problems in a sphere (e.g. the interior of a slowly rotating 
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Figure 4: The interplay between the kinetic energy -Ekin and magnetic energy 
Eb as well as the variation of the total energy during the creation of a toroidal 
component of the magnetic field from an initial purely poloidal magnetic field, 
due to differential rotation (see text for details). 



neutron star). The code works with spherical coordinates and spherical 
components of the velocity V and of the magnetic field B. The use of 
spherical coordinates and spherical components allows us to impose exact 
BCs on the surface of the star. Moreover, thanks to spectral methods, the 
singularities present on the z axis and at r — 0, which result from spher- 
ical coordinates, can be treated properly. We want to point out that the 
numerical solution is obtained in the whole star, while numerous dynami- 
cal spectral codes (for instance in geophysics) are restricted to some shells. 

As far as physics is concerned, order of magnitude considerations show 
that for most of the known compact stars, the acoustic time scale is much 
shorter than other timescales (Alfven, Brunt- Vaisala, rotation timescale, 
...). In order to easily have an efficient code, the anelastic approximation, 
which filters acoustic waves, is employed. In such a way, slow evolution of 
the star can be followed during tenths of its rotations. Moreover, due to 
the intrinsic accuracy and low dissipativity of spectral methods, the code 
is well suited to study the stability of steady-state MHD configurations 
and other physical problems that we shall investigate elsewhere. 
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